Precision shooting: Sampling long transition pathways 
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The kinetics of collective rearrangements in solution, such as protein folding and nanocrystal phase 
transitions, often involve free energy barriers that are both long and rough. Applying methods of 
transition path sampling to harvest simulated trajectories that exemplify such processes is typically 
made difficult by a very low acceptance rate for newly generated trajectories. We address this 
problem by introducing a new generation algorithm based on the linear short-time behavior of 
small disturbances in phase space. Using this "precision shooting" technique, arbitrarily small 
disturbances can be propagated in time, and any desired acceptance ratio of shooting moves can be 
obtained. We demonstrate the method for a simple but computationally problematic isomerization 
process in a dense liquid of soft spheres. We also discuss its applicability to barrier crossing events 
involving metastable intermediate states. 



I. INTRODUCTION 

Transition path sampling (TPS) is a versatile and ef- 
ficient set of computational techniques for the study of 
rare events J*2£& It has been successfully used to reveal 
the microscopic mechanisms of processes as diverse as au- 
toionization in liquid water—, structural transformations 
in nanocrystalline solids^, and folding of small proteins^ 
The purpose of this paper is to propose a new shoot- 
ing algorithm which can greatly increase the efficiency of 
TPS when transit times of activated trajectories greatly 
exceed the picosecond time scale of phase space stability. 

At its core TPS is a Monte Carlo procedure enabling 
a random walk in the ensemble of pathways that cross 
a free energy barrier between two metastable states (de- 
noted A and B). While this sampling is strongly biased 
towards reactive trajectories, it leaves the underlying dy- 
namics of the system unchanged. Thus, the result of a 
TPS simulation is a representative set of true dynami- 
cal pathways, weighted as if they were excerpted from 
an extremely long, unbiased simulation of equilibrium 
dynamics. Many analytical tools have been developed 
to extract from such a collection of trajectories useful 
molecular information about the process of interest. - 

The algorithm typically used to construct such a ran- 
dom walk is called shooting^. Here, a point along a 
given reactive trajectory is randomly selected and slightly 
changed; for instance, one might change the velocities 
of all particles by a small random number drawn from 
a symmetric distribution. Using the dynamical rules of 
the system, this shooting point is then propagated for- 
ward and backward in time to obtain a complete new 
trajectory. If this new trajectory still connects A with B 
it is accepted and used as a basis for the next shooting 
move; otherwise it is rejected. 

The efficiency of this algorithm in exploring the transi- 
tion path ensemble is based on a balance between the in- 



trinsic instability of complex dynamical systems and the 
local character of the shooting move: Small disturbances 
grow exponentially quickly in time, leading to separa- 
tion of trajectories typically within a few picoseconds. 
Nonetheless, if the disturbance is small, the new trajec- 
tory will be locally similar to the old one and is therefore 
likely to surmount the barrier between A and B; such 
shooting moves will be accepted frequently Just as with 
conventional Monte Carlo moves in configuration space, 
maximum efficiency can often be obtained by adjusting 
the size of the disturbance to achieve an acceptance prob- 
ability of roughly 40%^ 

Shooting moves are best suited for the study of sys- 
tems that relax quickly (within the picosecond time scale 
of trajectory separation) into their product state after 
reaching the top of the barrier. Many interesting pro- 
cesses, like the nucleation of first order phase transitions 
or conformational change in complex molecules, proceed 
much more slowly from the transition state. In TPS 
simulations of such systems, shooting moves must be 
made extraordinarily subtle in order to stand a reason- 
able chance of connecting reactant and product states. 
As a matter of practice, however, disturbances cannot be 
made arbitrarily small due to the limited machine pre- 
cision of floating point numbers. Lacking an ability to 
control the degree of global separation between trajecto- 
ries, TPS methods are severely compromised in efficacy. 
The demonstrated computational advantages of impor- 
tance sampling in trajectory space lose appeal when off- 
set by the wasted effort of generating a vast excess of 
non-reactive paths. 

In a recent paper-, Bolhuis addressed this problem 
by modifying slightly the rules that propagate a system 
in time. Specifically, a weak stochastic component was 
added to the dynamics, removing the unique correspon- 
dence between a trajectory's past and its future. It thus 
became possible to resample only parts of an existing 
pathway, leading to much higher acceptance probabili- 
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FIG. 1: Distance r, in reduced units, between two particles in 
a liquid of 108 soft spheres— as a function of time along a ref- 
erence trajectory and four shooting trajectories. At time zero, 
displacements of various size are added to particle velocities 
in the reference trajectory. Because phase space disturbances 
grow exponentially in time, decreasing the shooting displace- 
ment by successive orders of magnitude results in only linear 
increase in the time that elapses before trajectories separate. 
Shooting displacements smaller than 10 -15 can not be re- 
solved in double precision; the resulting shooting trajectories 
will exactly retrace the initial one. 



sequent time slices are related by 



X(i+l)At = $(»*At) 



(2) 



where the function $ propagates the system for one time 
step. 

Consider now a shooting move, in which a small dis- 
turbance Sxq is added to the shooting point xq to obtain 
state ya — xq + Sxq of the shooting trajectory Y . (To 
simplify notation, we will assume the shooting point to 
be Xq, the initial state of the trajectory, throughout this 
section. The algorithm we will describe applies trans- 
parently to shooting points at any chosen time along the 
trajectory.) Usually the perturbation 8x affects only 
momentum space, but changing the positions of the par- 
ticles can be useful in some cases^. The perturbed point 
2/o is then propagated for a number of n time steps to 
obtain y t — $ t (y ), where t = nAt and $ t refers to the 
n-fold application of the time step propagator. We de- 
fine the time-evolved disturbance Sxt by subtracting the 
old trajectory from the new one, Sxt = yt — %t- Due 
to the dynamic instability of the system, perturbations 
grow exponentially in time, 



ties for shooting moves and, Bolhuis reports, significant 
improvement in sampling efficiency^ 

In this paper, we show that it is possible to perform 
productive shooting moves for arbitrarily long transition 
paths without modifying a system's natural dynamics. 
Our technique for introducing and propagating extraor- 
dinarily small disturbances is based on the simple dynam- 
ics of small perturbations in phase space. We explain the 
method and offer a straightforward algorithm for imple- 
mentation in Sec. HT] Use of the technique is demon- 
strated in Sec. IIIII for a simple isomerization process in a 
dense liquid that, by construction, involves diffusive dy- 
namics on a rugged barrier. In Sec. IIVI we examine lim- 
itations of the method by considering reactive dynamics 
that pass through highly metastable, obligatory interme- 
diate states. 



II. LINEARIZED DYNAMICS OF SMALL 
PERTURBATIONS 

A. Exponential divergence of trajectories 

In a TPS simulation of a system evolving with deter- 
ministic dynamics, a trajectory X of length t consists of 
a number of "snapshots" XiAti which are separated by a 
time step At, 



\5x t \ « |&c |e 



XI 



(3) 



X = {x ,XAt,X 2 At, 



C T } 



(1) 



Here, the time slices XiAt are full phase space vectors, de- 
tailing the positions and velocities of all particles. Sub- 



Here, A is the largest Lyapunov exponent of the system*^ 
For typical fluid systems, 1/A « 1 ps. 

We wish to control precisely the time it takes for a 
small perturbation to reach a size of order 1, at which 
point the new trajectory will be essentially separated 
from the old one. This time determines the probabil- 
ity that the new trajectory will be reactive and therefore 
acceptable. Because of subsequent exponential growth, 
| foo| must be decreased by many orders of magnitude to 
increase the separation time of trajectories by even a few 
picoseconds (see Figs. [1] and [2]). With the standard dou- 
ble precision format for representing floating point num- 
bers on a computer, however, the smallest number that 
can be added to 1.0 to give a result distinguishable from 
1.0 is of the order of 10 -15 , and numerical results become 
unreliable at values of \Sxq\ well above this limit. (We 
assume throughout this paper that a system of units has 
been chosen such that typical numerical values of coordi- 
nates and momenta are of order 1.) Especially when the 
total length of the transition path is significantly longer 
than a few picoseconds, the limited range of practical 
displacement sizes constitutes a severe sampling prob- 
lem: Shooting moves will only be accepted from points 
in the vicinity of the barrier top; otherwise, new trajecto- 
ries will simply return to the stable state they came from 
and be rejected. As the system may stay near the a priori 
unknown barrier top only for a small fraction of the to- 
tal transition time, sampling can break down completely. 
In these cases, implementing shooting displacements of 
arbitrarily small size would be very helpful. 




Expanding y t around xq, we obtain 

5x t =Vt-x t = $t(*o + fan) - $t(x ) 



dx 



fa; + O [(8x ) 



(4) 



For small displacements \6xq\ < 10~ 15 , the linear approx- 
imation is for all practical purposes exact on the scale of 
a single time step, 



5x t = S 5xq 
where the matrix S is given by 



(5) 



FIG. 2: Time evolution of small displacements fat in a liq- 
uid of 108 soft spheres^ for disturbances of various sizes, 
\6xo\ — W~ a (legend values indicate values of a). All dis- 
placements grow exponentially with the same rate, up to 
the time where trajectories separate. The maximum possible 
value of \Sxt\ is determined by the dimensions of the simula- 
tion box. Note that adjacent lines are equidistant in the lin- 
ear regime, except for displacement sizes smaller than 1CP 12 . 
Although these smallest displacements yield trajectories that 
can in practice be distinguished from the base trajectory, lim- 
ited numerical precision introduces rounding errors that de- 
grade computational estimates of linear divergence. 
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FIG. 3: Time evolution of shooting moves from a point xo on 
the reference trajectory X. The two shooting displacements 
fao and fa , as vectors in high-dimensional phase space, 
point in the same direction but have different magnitudes, 
faj, = co fao- At a later time t short enough that first- 
order perturbation theory remains valid, these displacements 
remain proportional, 8x t = Co fat • The displacement of 
interest fat, no matter how small, can thus be constructed 
in the linear regime simply by dividing 8x t by the original 
scaling factor co. 



S = 



gjMgo) 

dxo 



(6) 



To integrate Sxq forward in time according to equation 
(O, the equations of motion for the matrix S could in 
principle be solved numerically. 12 Doing so in practice 
would be cumbersome, requiring calculation of all sec- 
ond derivatives of the interaction potential with respect 
to particle positions. We propose a much simpler ap- 
proach for advancing Sxt, inspired by methods for com- 
puting Lyapunov exponents in systems whose interaction 
potentials lack well-defined second derivatives. Our im- 
plementation is illustrated in Fig. [3] 

Instead of integrating the small perturbation Sxo, we 
follow the time evolution of a related perturbation Sxq , 
which is large enough to be added at the shooting point 
and propagated in the usual way, i.e., by integrating New- 
ton's equation of motion for y^ 



(i) _ 



X 



8xq ' . We use a 



superscript for j/W and 8x^ because in the following we 
will consider a family of different perturbed trajectories 



y(i 



- it® 



{y ,■■■,$■ }. with vl 



-,(i) 



Xi 



8x 



(i) 



Exploit- 



ing the linearity described by Eq. ([5|) , we choose 8£q ' to 
be in the same direction (in the high-dimensional phase 
space) as 8xq, 



8x q = co Sxq 



(7) 



:.(!> 



where ci is a scalar constant. If 8x is also small enough 
to justify the linear approximation of Eq. ([5]), 



SxP 



S Sx. 



(i) 



then the initial relationship between Sxq and Sx 
also at a later time t, 



(i) 



(8) 



holds 



B. Dynamics in the linear regime 



We propose to solve this problem by using perturbation 
theory to follow the time evolution of the displacement 
vector 8xq itself, up to the point where it grows large 
enough to allow accurate evaluation of the sum x t + 8x t . 



Sxt — S Sxq 



1 

CO 



SSx^ 



1 
co 



**i 1} 



(9) 



In the linear regime it is thus possible to follow the time 
evolution of arbitrarily small displacements Sxq by mon- 
itoring larger, proportional displacements. 

The linear approximation for the "helper" displace- 
in Eq. (JSJ) will of course remain valid for only 



ment 8x t 



*d 



a short time t\-' typically less than 1 ps. Our interest in 



the trajectory Y^\ however, is only as a proxy for the 
evolution of smaller displacements that cannot be repre- 
sented explicitly. As 5x t approaches the boundary of 
the linear regime, t ^ i lin , we may therefore switch our 
attention to a different helper trajectory Y^ 2 \ one whose 
displacement is initially too small to be of practical use 
but by the time t lin grows large enough to be represented 
explicitly. The new displacement 8x\ — y t — Xt can 
be obtained at any time t < £,. simply by scaling Sx t 



appropriately, 8x. 



(2) 



Sx t jc\. Because it is initially 



smaller than 8x\ , it will remain in the linear regime for 

a longer time, t\^ > t^. For times t > t\ h [ we there- 
fore proceed by integrating standard equations of motion 
for y( 2 ) until it approaches the boundary of the linear 
approximation. At that point we repeat the procedure, 
scaling back the displacement to switch attention to yet 
another helper trajectory. 

In effect we monitor a single displacement from the ref- 
erence trajectory whose magnitude is periodically scaled 
down such that the linear approximation is always valid. 
In this way we can monitor the time evolution of an ar- 
bitrarily small disturbance. The corresponding trajec- 
tory will be numerically indistinguishable from the refer- 
ence trajectory as long as the displacement's magnitude 
is smaller than ~ 10 -15 . At later times the displaced 
trajectory can be distinguished, and its dynamics can be 
safely computed by integrating equations of motion in 
the usual way. 




FIG. 4: Precision shooting algorithm for generating a trial 
trajectory Y (blue curve) whose initial displacement from the 
base trajectory X (red line) is extraordinarily small. Points in 
the light blue field (whose extent is « 10 _ ) cannot be nu- 
merically distinguished from the reference trajectory. Until 
the trial trajectory exits this region, its time evolution is cal- 
culated by proxy using "helper" trajectories Y^ 1 ' (thin black 
curves). Displacements of trial (black arrows) and helper tra- 
jectories from the base trajectory are related by proportion- 
ality as long as they remain within the linear regime, repre- 
sented by the light red field. To preserve this simple rela- 
tionship, the displacement of the helper trajectory Y^' from 
the reference trajectory is scaled back when it threatens to 
leave the region of linear dynamics, effectively switching the 
system to the next helper trajectory Y^' +1 ' . By following the 
systems dynamics along many sections of helper trajectories 
(green curve) and by keeping track of the rescaling factors, 
one can accurately construct the state of the trial trajectory 
once it becomes distinguishable from the reference trajectory. 
The result of this shooting move is a trial trajectory Y that is 
numerically identical to the base trajectory X over a certain 
length of time and then emerges from it in the correct way. 



C. Algorithm 

This insight suggests the following algorithm, which 
implements a shooting trajectory Y, whose initial devia- 
tion Sxq from the base-trajectory X is smaller than the 
precision limit. This is done by monitoring "helper" tra- 
jectories jw') , which are obtained by repeated rescaling. 
For an illustration of this algorithm see Fig. 2] 

1. At the shooting point xq, add a displacement 5x 

of fixed size \8x \ = a. The displacement 8x Q is 
parallel to Sxo and larger by a factor of cq. 

2. Propagate the point j/q = xq + 8x forward in 
time for n time steps, corresponding to a time in- 
terval t = nAt. 

3. Compute the factor c\ — \8x\ \j\8xq \ — \y\ '. — 

Xt\/\8x | quantifying the divergence from the ref- 
erence trajectory. Switch to a new helper trajectory 
by setting y\ = Xt + 8x t /ci = x t + 8x t . Store 
the factor c\. 

4. Propagate the new displacement forward in time 
by integrating the equations of motion for n steps 
beginning from y t . Calculate and store the factor 



-,<J) 



and 

At 



C2 



\8x%\l\8xf\ 



5. Iterate step 4, each time beginning from y, ., 

rescaling by the factor Cj = \8x , 3 t \/\8x, 3 _ x ^ t \. 
every iteration compute the displacement of inter- 
est, Sxjt — Cj 1 5Xj t , where Cj = n|,~ Cfc is the 
product of all factors used for rescaling so far. To 
store the current point along the actual shooting 
trajectory, compute y Jt = Xjt + Sxjt- As long as 
\$ x jt\ ^ 10~ 15 , yjt will be numerically identical to 
x jt . 

6. If \Sxjt\ > c, the actual shooting displacement is 
large enough to be treated in the usual way: Set 
yjt — Xjt + Sxjt and integrate equations of motion 
from this point without further rescalings (ceasing 
iteration of step 4). 

Note that for shooting moves conducted at points other 
than xq, the procedure must be repeated backward in 
time to obtain a complete shooting trajectory. In the 
following we discuss the accuracy of this scheme and give 
recommendations for choosing values of a and n. 



D. Validity of the linear approximation 

The above algorithm is exact only if the linear ap- 
proximation of Eq. ([5]) holds. For perturbations of fi- 
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FIG. 5: Time evolution of the magnitude of two shooting 
displacements for a fluid of WCA particles.— Displacements 



of size \Sx { L '\ = 1CT 7 (black) and \8x ( ">\ = 10" 6 (red) are 
initially proportional (pointing in the same direction in phase 
space). Both are rescaled to their initial size every 100 time 
steps (see inset for a magnified view) to preserve this linear 
relationship. Deviation from proportionality is quantified by 
the relative error e(t) (blue) defined in Eq. (1101 ) . 



nite size, deviations from this approximation occur. The 
question thus arises, how accurate an approximation is 
this approach for propagating small disturbances? More 
specifically, to what extent do helper displacements re- 
main proportional to the actual shooting displacements 
of interest? One could certainly imagine that the fast 
growth of small non-linearities rapidly erodes the linear 
relationship on which we depend. Here we present evi- 
dence from computer simulations that proportionality of 
small displacements can hold in practice over very long 
time scales. 

Figure [5] shows the time evolution of two proportional 
disturbances. The initial displacement vectors point in 
the same direction of phase space but have different mag- 



nitudes, \5£q J 



= 1(T 7 and \6x% 



(2) 



10" 



Thc 



spective shooting trajectories were propagated indepen- 
dently, and the displacements from the base trajectory 
were rescaled to their initial length every 100 time steps. 
For a perfectly linear time evolution, these displacements 
remain proportional at all later times. In practice the 
vectors Sx). and 5x t will develop a nonzero angle due 
to non-linearities. To quantify this deviation from paral- 
lel alignment, we define 



e(t) 



|&rf ] - SxP/lO] 



Sx 



(i) 



(10) 



As shown in Fig. [5l the relative error e(t) does not grow 
above a low level even for very long simulation runs. 

This long time stability of aligned disturbances holds 
over a broad range of displacement sizes between 10~ 10 
and 10 -3 . Values of e(t) can be somewhat larger than for 
the specific case plotted in Fig. \E\ but on average do not 
grow larger than 10 -3 for any case. The error is insensi- 
tive to the choice of the rescaling interval nAt, as long as 
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FIG. 6: Relative orientations of three phase space displace- 
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ment vectors (Sx). 1 '' , 8x) z> , and 8x\ ) for a WCA fluid.— 

We plot the angle 7ij - = cos _1 (W| • ^Vl^II^Fl) be- 
tween each pair of displacements as a function of time. Dis- 
placements differ in initial direction and size (\8x Q | = 10 -8 , 
10~ 8 , \8x Q | = 10~ 7 ) and are rescaled to their origi- 



i*4 2)| 



nal sizes at different intervals (n x — 100, n y — 200, n z — 500). 



the displacements stay smaller than approximately 10~ 2 , 
where the linear regime breaks down. For displacements 
smaller than 10 -10 , rounding errors become problematic, 
and phase space displacements do not remain parallel to 
a good approximation. One might expect the region of 
long time stability to extend to even lower levels, closer 
to the precision limit of 10~ 15 . However, the total achiev- 
able accuracy in a computer simulation depends, among 
other factors, on the details of the integrator $, the size 
of the time step, and the dimensionality of the system, 
and can lie well above the precision limit of 10~ 15 . We 
find that for the particular system used here, the preci- 
sion level is about 10~ 12 . 

In the light of these observations, a value of a = 10~ 6 
as initial magnitude for helper displacements seems ap- 
propriate. This choice lies midway between the upper 
limit of the linear regime (approximately 10~ 2 ), and the 
point where rounding errors become dominant (approxi- 
mately 10~ 10 ). As the accuracy of the rescaling scheme 
is quite insensitive to the frequency of rescalings, many 
equally good choices of n are possible. As a starting 
point, a value of n which leads to rescalings every time 
the displacements have doubled their size is advisable. 

The fact that e(i) does not show any systematic long- 
time growth in Fig. [5] seems surprising. After all, no 
constraint is imposed on the direction of the displace- 
ment vectors. Why does an accumulation of errors not 
eventually lead to decoupling and e(t) s=s 1? Stability 
of the precision shooting algorithm is in fact a simple 
and direct consequence of the collective dynamics of dis- 
placements in the linear regime. In Figure [6] we plot 
the angles between three periodically rescaled shooting 
displacement vectors of different size and random initial 
direction. Eventually, they all rotate into the same di- 
rection, which is associated with the largest Lyapunov 
exponent A of the system. The time scale on which the 



directions of different displacement vectors converge is on 
the order of 1/AA, where AA is the difference between 
the first and second largest Lyapunov exponents.— It is 
because of this convergence, that the difference vector 
Sx t — 8x\ ' jc between two proportional displacements 
with initially identical direction will stay small. 

We point out that this property constitutes a main 
difference of our method over the stochastic scheme in- 
troduced by Bolhuis- and similar algorithms. Consider, 
for instance, the following simple algorithm that can be 
viewed as a smooth version of the stochastic scheme by 
Bolhuis: 

• Choose a shooting point a; s At- 

• A fixed number of timesteps n earlier and later, at 
the points X( s +„)At an d Z( s -n)Ati add a displace- 



ment of 10" 
particle. 



-is 



to one velocity component of one 



• Integrate the points £( s + n )At and a;( s _ n )A4 forward 
and backward in time, respectively, to get a com- 
plete new trajectory. 

Just like the precision shooting algorithm, this simple 
scheme results in a shooting trajectory that is numeri- 
cally identical to its base trajectory for a certain period 
of time. However, the emerging separation between base 
and shooting trajectories will not be consistent with a 
shooting move conducted at xq, but rather with two un- 



corrected shooting moves at x 



(s+n)At 



and . 



: (s-n)At- 



Oi 



algorithm, on the other hand, correctly reproduces the 
correlated forward and backward dynamics of a displace- 
ment introduced at x s At- 



III. A SIMPLE TEST SYSTEM 

We demonstrate the precision shooting algorithm on 
a simple isomerization process of a solvated diatomic 
molecule in three dimensions. 

Our test system consists of 389 particles interacting 
via the WCA-i^ potential. We use conventional reduced 
units, with particle mass and potential parameters a and 
e all set to unity. Particles #1 and #2 do not interact 
via the WCA potential, but are bonded through a one- 
dimensional potential with two deep minima separated 
by a rough barrier (see Fig. [TJ: 



v(x) = 



Here, q(x) 




- q(x) 2 /w 2 ] 

-{q{x)-b) 2 /w 2 } 2 
h cos 2 [a(q(x)-b/ 2)} 
2 y/l+ga 2 ( q (x)-b/2)2 



if«(aO<0, 

if q(x) > b , 
else. 

(11) 
(r c + w), x = X2 — x\ is the difference 
between the a;-component of the position of the bonded 
particles, r c — 2 1 ' 6 is the cutoff of the WCA potential, 
w = 1 determines the width of the minima, b = 10 and 
hi = 10 are the length and height of the barrier in be- 
tween, respectively, and the constants hi = 3, a = 77r/6, 




FIG. 7: Potential energy v(x) of interaction between particles 
#1 and #2, comprising the diatomic molecule in our model 
isomerization process, plotted as a function of the difference 
x = X2 — xi between their ^-coordinates. The dashed lines 
mark the boundaries of the minima A and B. 
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FIG. 8: Distribution of transition times T for our model iso- 
merization process, as gauged from 2500 trajectories initiated 
at the barrier top. Inset: Difference a; = X2 — xi between the 
^-coordinates of particles #1 and #2 as a function of time for 
a typical trajectory. 



and 5 = 2 determine the shape of the barrier. The po- 
tential and its derivative are continuous by construction. 

To speed computation, we borrow a trick from Bolhuis' 
work: 8 Particle #2 is considered to lie always to the right 
of particle #1, hence x > 0. This choice, together with 
the one-dimensionality of v (x), allows us to choose a sim- 
ulation box with dimensions 14.4 x 6 x 6. The resulting 
particle density is 0.75, the total energy per particle is 
1.0, and the temperature is 0.45, as gauged by average 
kinetic energy. We use the velocity Verlet algorithm 14 
to integrate the equations of motion with a time step of 
0.002. 

We are interested in sampling the transition of the 
dimer from the "contracted" minimum A at xa = r c 
to the "extended" minimum B at xb = r c + b + 2w. The 
dimer is defined to be in state A for x < xa + 0-75w 
and in state B for x > xb — 0.75u> (see Fig. [7]). Be- 
cause the system is dense, and the barrier is both long 
and rough, relaxation from the transition state into either 



stable minimum is quite protracted. 

In conducting TPS simulations it is important that 
sampled trajectories are not shorter than typical sponta- 
neous barrier-crossing events^ We determine this typical 
duration for our simple model system by initiating many 
straightforward molecular dynamics simulations with the 
dimer bond length set at x = r c + w + b/2, corresponding 
to the middle of the barrier. Integrating the equations 
of motion forward and backward in time yields a repre- 
sentative sample of the transition path ensemble. For a 
particular trajectory, the transition time T is the time 
the system spends between regions A and B. The result- 
ing distribution of transition times is plotted in Fig. [8] 
For TPS simulations, we choose a total trajectory length 
of 3 x 10 5 time steps, long enough to include 98% of the 
natural transition path ensemble. The bias of our sam- 
pling to short transitions is therefore minor. 

Although the artificial potential energy landscape 
studied here does not directly represent any physical sys- 
tem of interest, it nevertheless shares with many real 
systems features that lead to long transition pathways 
and make straightforward application of TPS methods 
ineffective. In our view roughness of the barrier region 
is an important ingredient. Models featuring long but 
flat barriers, such as that of Ref4, should not in fact 
pose any severe problems for path sampling via the stan- 
dard shooting move. Assuming that motion atop such a 
flat barrier is diffusive in nature, and that time evolution 
from the edge of the barrier proceeds into the adjacent 
minimum with near certainty, then a trajectory initiated 
on the barrier will relax into stable state A with proba- 
bility pa = 1 — y/b, where y is the initial distance from A 
and b is the width of the barrier. Similarly, the probabil- 
ity of relaxing first into state B is pe = y/b. A standard 
shooting move from the barrier region then yields a re- 
active trajectory with probability 
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This value of the acceptance rate should correspond to 
near optimal sampling of the transition path ensemble. 2 
A problematically low acceptance rate would only arise 
if one were to sample trajectories of insufficient length, 
i.e., paths shorter than typical spontaneous transitions. 
In our TPS simulations, only momenta (and not parti- 
cle positions) are disturbed in the shooting moves, with 
each particle's momentum changed in each direction by 
an amount drawn from a Gaussian distribution of stan- 
dard deviation Ap (followed by rescaling of all momenta 
to enforce energy conservation) X We conduct standard 
shooting moves with values of Ap of 1CP 1 , 10~ 5 , and 
10~ 10 , as well as precision shooting moves with Ap rang- 
ing in size from 1CP 10 to 10~ 300 . The latter are im- 
plemented using helper displacements with Ap = 1CP 7 
and are rescaled every time they reach twice their orig- 
inal size. For the system size studied here, the ini- 
tial magnitudes of the resulting displacement vectors are 
larger than the corresponding value of Ap by a factor 
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FIG. 9: Fraction of shooting moves for our model isomer- 
ization process that are accepted in long TPS simulations. 
Acceptance ratios are shown for shooting displacements of 
various sizes, Ap = 10 _a , implemented using the standard 
shooting algorithm for values of a of 1, 5, and 10, and the 
precision shooting algorithm for a > 10. For a — 10, the 
result obtained from the precision shooting algorithm (red) is 
effectively indistinguishable from the one obtained with stan- 
dard shooting. 



of roughly 34 on average. For instance, Ap = 10 -7 re- 
sults in displacement vectors with an initial size of about 
| fool ~ 3.4 x 10~ 6 . For each set of sampling parame- 
ters, we attempt 50,000 Monte Carlo moves in trajectory 
space. Roughly half of these trial moves are generated by 
shooting. The other half are generated by a procedure 
called "shifting" j 2 - in which short trajectory segments are 
added to and subtracted from the ends of an existing 
path. 

Figure [9] shows the fraction of attempted shooting 
moves that are accepted in TPS simulations of the di- 
atomic isomerization with a rough barrier. While stan- 
dard shooting moves are accepted with low frequency, 
any desired acceptance ratio can be obtained by using the 
precision shooting technique. Figure [TTJ] shows changes in 
transition time T over the course of two TPS runs with 
shooting displacements of Ap = 10 _1 and Ap = 10~ 100 . 
A dramatic difference in the efficiency of generating qual- 
itatively different trajectories for the two cases is evident. 

To assess the improvement in sampling efficiency 
achieved with precision shooting, we quantify the compu- 
tational effort necessary to generate statistically indepen- 
dent transition pathways. More specifically, we calculate 
the autocorrelation function 
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where ST(n) — T(n) — (T) is the deviation of the tran- 
sition time after the n-th shooting move from its aver- 
age (T), as calculated from all collected trajectories.^ 
Rapid decay of c(n) indicates an efficient sampling of 
trajectories. Figure [TT] shows the logarithm of c(n) for 
different shooting displacement magnitudes along with 
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FIG. 10: Variation in isomerization transition time T over the 
course of long TPS runs. The bold black line shows results for 
a simulation using shooting displacements with Ap = 10 - , 
while the thin green line corresponds to Ap = 10~ 100 . 
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FIG. 11: Number of shooting moves required to generate a 
statistically independent isomerization trajectory using shoot- 
ing displacements of various sizes, with Ap = 10~ a . Inset 
shows decay of correlation c(n) in the transition time T follow- 
ing n attempted shooting moves. The "decorrelation time" v 
is defined as the value of n beyond which c(n) < 1/2. 



the "decorrelation time" v, denned as the number of suc- 
cessive shooting moves after which the correlation func- 
tion decays to a value less than 1/2. The maximal sam- 
pling efficiency is achieved for shooting displacements 
with Ap i=s 10~ 100 . Improvement over the largest dis- 
placement we considered (Ap = 10 _1 ) is more than ten- 
fold. Following Bolhuis 8 , we also investigate as a measure 
of decorrelation changes in the bond length x midway in 
time through the crossing event. Decay of correlations 
in this quantity, and the implied dependence of sampling 
efficiency on shooting displacement size, mirror those re- 
ported for the transition time T . 

As Fig. [11] illustrates, sampling is comparably effi- 
cient for a broad range of displacement sizes between 
Ap = 1(T 60 and Ap = 1CT 260 . In this regime, the effi- 
ciency gain due to increased acceptance rates for smaller 
shooting moves is compensated almost exactly by the effi- 
ciency loss due to increased similarity between the shoot- 
ing trajectory and its base trajectory. Using equation 
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FIG. 12: TPS simulation of isomerization dynamics which 
must proceed through a deep intermediate energy minimum. 
Top left: Interaction potential between particles #1 and #2. 
Top right: Acceptance ratio of shooting moves as a function 
of displacement size, Ap = 10 _O! . Bottom: Bond length x 
as a function of time for the starting trajectory and for a 
pathway obtained after many shooting and shifting moves. 
The two trajectories have been shifted in time to highlight 
their similarity in the vicinity of the intermediate state. 



([3]), the time T^ over which a shooting trajectory with 
displacement size 10~ Q cannot be resolved from its base 
trajectory can be approximated by the time required for 
the displacement size to reach 10 ~ 15 , 
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For our system, 1/Ai w 150 Ai and therefore T id rs 10 5 
time steps for the smallest displacement with Ap = 
10 -300 . Even for this small displacement size, Xld is 
only 30% of the total trajectory length L and efficient 
sampling is still possible. If Sxq is decreased further, 
Tjd will become comparable to L and sampling effi- 
ciency will decrease accordingly. 1 For a displacement 
with Ap = 10~ 60 , the largest size that leads to optimum 
efficiency, T|d amounts to only 5% of the total trajectory 
length. 



IV. METASTABLE INTERMEDIATE STATES 

By extending the time span over which a shooting tra- 
jectory tracks its base trajectory, the algorithm proposed 
in this work can substantially increase the efficiency of 
TPS simulations that suffer from poor acceptance of 
shooting moves. The method is fully consistent with de- 
terministic dynamics and faithfully reproduces the diver- 
gent behavior of arbitrarily small displacements in phase 
space. We emphasize, however, that the method does 
not solve all problems whose primary symptom is a low 
shooting acceptance rate. Most importantly, it does not 
overcome challenges associated with metastable interme- 
diate states. In this section we explore the this difficulty 
in the context of diatomic isomerization. 



In order to explore the consequences of metastable 
intermediates, we have modified the diatomic potential 
v(x) to include a deep minimum midway between con- 
tracted and extended states (see Fig. [l~2")) . 

f h^l-qixf/wl] 2 if g(x) < 0, 

«(«) = < hi [1 - (?(ar) - 2 W2 ) 2 / W 2] 2 ^ if g(x) > 2w2j 
[ ft-i - /i 2 [l - (q(x) - w 2 ) 2 /wl] 2 else. 

(15) 
Here, q(x) — i-(r c +wi), w\ — 0.5, 102 = 3, h\ = 15, and 
h'2 = 10. Limited by machine precision, standard shoot- 
ing moves fail completely in this case: Even shooting 
moves initiated near the intermediate minimum C rapidly 
separate from their base trajectories and with high prob- 
ability do not escape to stable state A or B. Only with 
the precision shooting technique, using a displacement 
size smaller than 10~ 20 , are we able to conduct successful 
shooting moves. This success does not indicate, however, 
that trajectory space is sampled efficiently: A compari- 
son of the first trajectory^ with a pathway obtained af- 
ter many thousands of shooting moves shows that those 
parts of the trajectory spent within the intermediate are 
not resampled at all (see Fig. |T2j) ; they are numerically 
identical. 

Transitions involving strongly metastable intermedi- 
ates are in fact fundamentally problematic for TPS meth- 
ods, unless the dynamics of intermediates' appearance 
and disappearance can be identified as distinct kinetic 
substeps. If the typical time spent in C is manageable in 
a computer simulation, then the intermediate does not 



pose a problem even to the standard shooting move. If, 
on the other hand, the free energy barriers delimiting the 
intermediate state are large compared to typical thermal 
excitations, then escaping C will itself be a rare event. 
In such cases, typical transitions from A to B require 
at least two unlikely fluctuations (activating entry and 
exit of each intermediate state), well separated in time. 
Any shooting move that perceptibly modifies dynamics 
between these rare fluctuations will be rejected with high 
probability. Precision shooting can readily generate sub- 
tly modified pathways that remain reactive but cannot be 
expected to effectually switch between reactive trajecto- 
ries that follow substantially different courses through the 
intermediate state. As TPS leaves a system's natural dy- 
namics unchanged, it can eliminate only the largest time 
scale associated with a rare event. Without resorting to 
methods that prescribe in some sense the detailed route 
between stable states, one can overcome the challenge 
of metastable intermediates with TPS only by subdivid- 
ing transition dynamics into several steps, each of which 
involves a single dynamical bottleneck. 
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